;********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;load library for  wetbulb_stull
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/heat_stress.ncl"
;load "./panel_two_sets.ncl" ; have to download this script to run
;********************************************
 begin

;set baseline
 dirPath = "~/source_data/"
 dirPath2 = dirPath+"/outdoor/"

 mo  = (/"ACCESS-CM2","AWI-CM-1-1-MR","BCC-CSM2-MR","CMCC-CM2-SR5","EC-Earth3","GFDL-ESM4", \
            "IITM-ESM", "KIOST-ESM","MIROC6","MIROC-ES2L","MPI-ESM1-2-HR","MRI-ESM2-0"/);

 case = (/"ssp585" /)
 ncase = dimsizes(case);
 nmod = dimsizes(mo)
 levs = (/1,1.5,2,2.5,3,3.5,4/); 7 warming levels
 nlev = dimsizes(levs)

 syr = 1990
 fyr = 2099;
 csyr = sprinti("%0.4i",syr);
 cfyr = sprinti("%0.4i",fyr);
 nyr = fyr - syr + 1;
 ndays = nyr*365

 ;read dimension for 1 deg data
  dimPath=dirPath+"/domain/"
  landfrac = "sftlf_fx_360x180_1deg_gn.nc"
  area = "areacella_fx_360x180_1deg_gn.nc"

  fdim  = addfile(dimPath+landfrac, "r")
  fdim2 = addfile(dimPath+area, "r")
  lat := fdim->lat
  lon := fdim->lon
  nlat = dimsizes(lat);
  nlon = dimsizes(lon);
  garea = fdim2->areacella  ;m^2
  lfrac = fdim->sftlf/100; %
  wgt = garea*1e-6*lfrac ; sqrkm
  wgt@_FillValue = lfrac@_FillValue; make sure FillValue is set correctly for wgt
  copy_VarCoords(garea, wgt)
  printVarSummary(wgt)
  print("global total land area: "+sum(wgt)*1e-6+"million sq km")



;===Use precalculated global warming amount relative to 1980-2009 (based on ERA5, HadCRUT5, BerkeleyEarth, NOAAGlobalTemp) to classify warming levels for each CMIP6 model
  ;see Data_global-warming-amount-ERA5-30yr-CMIP6_12models.ncl
  filename1 = dirPath+"Global_warming_amount_10yrma_"+syr+fyr+"_12models_obs30yr.csv"; 10-year running mean
  filename2 = dirPath+"Global_warming_amount_30yrma_"+syr+fyr+"_12models_obs30yr.csv"; 20-year running mean

;---Read in file as array of strings so we can parse each line
  lines  = asciiread(filename1,-1,"string")
  lines2 = asciiread(filename2,-1,"string")
  k = 1 ;number of header lines
  nlines = dimsizes(lines)-k   ; skip header lines
  delim = ","
;---Read fields (columns)
  models  =          str_get_field(lines(k:),1,delim)
  cases   =          str_get_field(lines(k:),2,delim)
  central_10yr =  tointeger(str_get_field(lines(k:),3,delim))
  period10yr =       str_get_field(lines(k:),5,delim)
 ;read mean temperature anomaly from file 1
  Ta_anomaly10yr = tofloat(str_get_field(lines(k:),7,delim))
  Ta_anomaly10yr@_FillValue=1e+20

  central_30yr =    str_get_field(lines2(k:),3,delim)
  period30yr =       str_get_field(lines2(k:),5,delim)
  Ta_anomaly30yr = tofloat(str_get_field(lines2(k:),7,delim))
  Ta_anomaly30yr@_FillValue=1e+20


 ;------choose whether 10-yr or 30-yr running period
  Ta_anomaly = Ta_anomaly30yr; either 10-yr or 30-yr running median/max
  periods = period30yr
  central_year = central_30yr

  levels := fspan(1,4,7); 1 to 4 degC by 0.5 degC
  nlev = dimsizes(levels)


;==constant model parameters==
 alpha=0.3  ;reflectance of human body depending on clothes
 rho=1.2 ;density of air (kg/m3)
 cp=1010 ;specific heat capacity of air (Joules/kg/°C)
 ;=baseline skin temperature
 Ts=36 ;skin temperature in hot conditions, beyond which core temperature will rise (see Supplementary Table 1)
 esat_skin = satvpr_water_bolton(Ts, (/0,1/))  ;saturation vapor pressure at skin temperature 
 ;=radiation terms
 e_skin = 0.97 ;skin emissivity (Brewster, M. Quinn (1992))
 fs = 0.8  ;the fraction of total skin area effectively involved in radiant heat exchange (Steadman 1979, Fiala et al. 1999)

 ;=metabolic heat
 ;mean human skin surface area 1.7 (1.6-1.8 m^2; Parsons 1999)
 ;total metabolic heat generation 300 W for moderate work (Kjellstrom 2013)
 MH = 300/1.7  ;outdoor moderate work (W m^-2)



 ;***below module calculates the diel cycle of G at local time of each grid cell using concurrent climate variables when annual maximum Gday happens
 
 DATA = True; 
 if (DATA) then
 
  nstep = 8; steps per day; 3-hrly data
  
  Gday_sun_ens = new((/nmod,nlev,nlat,nlon/),"float")
  Gday_diff_ens = new((/nmod,nlev,nlat,nlon/),"float")
  Gday_zero_ens = new((/nmod,nlev,nlat,nlon/),"float")
  G_sun_ens = new((/nmod,nlev,nstep,nlat,nlon/),"float")
  G_diff_ens = new((/nmod,nlev,nstep,nlat,nlon/),"float")
  G_zero_ens = new((/nmod,nlev,nstep,nlat,nlon/),"float")
  Z_sun_ens = new((/nmod,nlev,nstep,nlat,nlon/),"float")
  Z_diff_ens = new((/nmod,nlev,nstep,nlat,nlon/),"float")
  Z_zero_ens = new((/nmod,nlev,nstep,nlat,nlon/),"float")
  do mm=0,nmod-1

   print(" ")
   print("process CMIP6 model: "+mo(mm))


    ;==scenario 1: full solar radiation outdoors
    f0 = addfile(dirPath2+"/G-daynight-"+syr+fyr+"-"+mo(mm)+"-vcbc-outdoor.nc", "r")
    Gday_sun  := f0->Gday_mx_sun;
    Gday_diff  := f0->Gday_mx_diff;
    Gday_zero  := f0->Gday_mx_zero;

   nc1 = dirPath2+"/Gday_G-Z-3h-"+syr+fyr+"-"+mo(mm)+"-vcbc-outdoor.nc"
   if (fileexists(nc1)) then
    print("read 3hrly G and sun zenith angle")
    f1 = addfile(nc1, "r")
    G_sun  := f1->G_Gday_sun
    Z_sun  := f1->Z_Gday_sun    
    G_diff  := f1->G_Gday_diff
    Z_diff  := f1->Z_Gday_diff    
    G_zero  := f1->G_Gday_zero
    Z_zero  := f1->Z_Gday_zero 
   
   else
    f1 = addfile(dirPath2+"/Gday_sun-TQRUPZ-"+syr+fyr+"-"+mo(mm)+"-vcbc-outdoor.nc", "r")
    taC := f1->Ta_Gday_sun; (yrs:yrd,:,:) ; degC;
    rh  := f1->RH_Gday_sun; (yrs:yrd,:,:) ; relative humidity
    pa  := f1->P_Gday_sun;
    Rs  := f1->Rs_Gday_sun; (yrs:yrd,:,:); total solar radiation (daily)
    Ld  := f1->Ld_Gday_sun; (yrs:yrd,:,:); downward longwave at surface (daily)
    Lg  := f1->Lg_Gday_sun; (yrs:yrd,:,:);  upwelling longwave from ground
    Rd  := f1->Rd_Gday_sun; (yrs:yrd,:,:);  diffuse radiation
    Rg := f1->Rg_Gday_sun; (yrs:yrd,:,:); reflected shortwave from ground
    wind2m := f1->U2_Gday_sun; (yrs:yrd,:,:); 
    Rsa  := f1->Rsa_Gday_sun;

    ;===use Fiala convective heat transfer function
    hc:= 14.1*wind2m^0.5 ;forced convection function from Fig.3 in Fiala 1999
    hc:= where(wind2m .le. 0.1, 3.3, hc); set natural convection hc=3.3 W m-2 K-1 for U<=0.1 m/s as in Fiala

    ;==saturation vapour pressure==
    esat_skin = satvpr_water_bolton(Ts, (/0,1/))  ;esat at skin temperature, degC
    esat_air  := satvpr_water_bolton(taC, (/0,1/)) ;esat at air temperature, degC to Pa
    e_a := (rh/100)*esat_air
   
    ;=longwave radiation input (from air and ground) and output (emitted by skin)
    RLin := e_skin*0.5*(Ld + Lg)
    RLout = e_skin*5.67*10^-8*(Ts+273.15)^4 ;longwave from skin

    ;Latent heat of vaporization as a function of temperature==
    ;should use skin and sweat temperature, not air temperature
    lambda = latent_heat_water(Ts, (/0, 0/), 1, False)
    gamma0 = (pa*cp)/(0.622*lambda)

    H  := hc*(Ts - taC)
    LE := hc/gamma0*(esat_skin - e_a)

    LEmax = 500 ; [W m-2] or 1250 g/h for acclimatized persons according to ISO 7933:2023
    ;LEmax = 400 ; [W m-2] or 1000 g/h for unacclimatized persons according to ISO 7933:2023
    LE := where(LE .gt. LEmax, LEmax, LE);
 
    Rn := RLin - RLout + Rsa
    ;add an effective radiation area factor (Steadman 1979a)
    Rn = fs*Rn ; apply the effective radiant area factor fs=0.8 (Steadman 1979a) 

    G_sun := Rn - H - LE + MH
    copy_VarCoords(taC, G_sun)


    ;=scenario 2: exposure only to unavoidable diffuse radiation (under shading objects)

    f1 = addfile(dirPath2+"/Gday_diff-TQRUPZ-"+syr+fyr+"-"+mo(mm)+"-vcbc-outdoor.nc", "r")
    taC := f1->Ta_Gday_diff; (yrs:yrd,:,:) ; degK;
    rh  := f1->RH_Gday_diff; (yrs:yrd,:,:) ; relative humidity
    pa  := f1->P_Gday_diff;
    Ld  := f1->Ld_Gday_diff; (yrs:yrd,:,:); downward longwave at surface (daily)
    Lg  := f1->Lg_Gday_diff; (yrs:yrd,:,:);  upwelling longwave from ground
    Rd  := f1->Rd_Gday_diff; (yrs:yrd,:,:);  diffuse radiation
    Rg := f1->Rg_Gday_diff; (yrs:yrd,:,:); reflected shortwave from ground
    wind2m := f1->U2_Gday_diff; (yrs:yrd,:,:);
    Z_diff := f1->Z_Gday_diff;

    ;===use Fiala convective heat transfer function
    hc:= 14.1*wind2m^0.5 ;forced convection function from Fig.3 in Fiala 1999
    hc:= where(wind2m .le. 0.1, 3.3, hc); set natural convection hc=3.3 W m-2 K-1 for U<=0.1 m/s as in Fiala

    ;==saturation vapour pressure==
    esat_skin = satvpr_water_bolton(Ts, (/0,1/))  ;esat at skin temperature, degC
    esat_air  := satvpr_water_bolton(taC, (/0,1/)) ;degC to Pa
    e_a := (rh/100)*esat_air

    ;=longwave radiation input (from air and ground) and output (emitted by skin)
    RLin := e_skin*0.5*(Ld + Lg)
    RLout = e_skin*5.67*10^-8*(Ts+273.15)^4 ;longwave from skin

    ;Latent heat of vaporization as a function of temperature==
    ;should use skin and sweat temperature, not air temperature
    lambda = latent_heat_water(Ts, (/0, 0/), 1, False)
    gamma0 = (pa*cp)/(0.622*lambda)

    H  := hc*(Ts - taC)
    LE := hc/gamma0*(esat_skin - e_a)

    LEmax = 500 ; [W m-2] or 1250 g/h for acclimatized persons according to ISO 7933:2023
    LE := where(LE .gt. LEmax, LEmax, LE);

    Rn2 := RLin - RLout + 0.5*(1-alpha)*(Rd+Rg)
    ;add an effective radiation area factor (Steadman 1979a)
    Rn2 = fs*Rn2 ; apply the effective radiant area factor fs=0.8 (Steadman 1979a) 

    G_diff := Rn2 - H - LE + MH
    copy_VarCoords(taC, G_diff)

    ;=scenario 3: exposure only to longwave radiation (completely dark)

    f1 = addfile(dirPath2+"/Gday_zero-TQRUPZ-"+syr+fyr+"-"+mo(mm)+"-vcbc-outdoor.nc", "r")
    taC := f1->Ta_Gday_zero; (yrs:yrd,:,:) ; degK;
    rh  := f1->RH_Gday_zero; (yrs:yrd,:,:) ; relative humidity
    pa  := f1->P_Gday_zero;
    Ld  := f1->Ld_Gday_zero; (yrs:yrd,:,:); downward longwave at surface (daily)
    Lg  := f1->Lg_Gday_zero; (yrs:yrd,:,:);  upwelling longwave from ground
    wind2m := f1->U2_Gday_zero; (yrs:yrd,:,:);
    Z_zero := f1->Z_Gday_zero;

    ;===use Fiala convective heat transfer function
    hc:= 14.1*wind2m^0.5 ;forced convection function from Fig.3 in Fiala 1999
    hc:= where(wind2m .le. 0.1, 3.3, hc); set natural convection hc=3.3 W m-2 K-1 for U<=0.1 m/s as in Fiala

    ;==saturation vapour pressure==
    esat_skin = satvpr_water_bolton(Ts, (/0,1/))  ;esat at skin temperature, degC
    esat_air  := satvpr_water_bolton(taC, (/0,1/)) ;degC to Pa
    e_a := (rh/100)*esat_air

    ;=longwave radiation input (from air and ground) and output (emitted by skin)
    RLin := e_skin*0.5*(Ld + Lg)
    RLout = e_skin*5.67*10^-8*(Ts+273.15)^4 ;longwave from skin

    ;Latent heat of vaporization as a function of temperature==
    ;should use skin and sweat temperature, not air temperature
    lambda = latent_heat_water(Ts, (/0, 0/), 1, False)
    gamma0 = (pa*cp)/(0.622*lambda)

    H  := hc*(Ts - taC)
    LE := hc/gamma0*(esat_skin - e_a)

    LEmax = 500 ; [W m-2] or 1250 g/h for acclimatized persons according to ISO 7933:2023
    LE := where(LE .gt. LEmax, LEmax, LE);

    Rn3 := RLin - RLout
    ;add an effective radiation area factor (Steadman 1979a)
    Rn3 = fs*Rn3 ; apply the effective radiant area factor fs=0.8 (Steadman 1979a) 

    G_zero := Rn3 - H - LE + MH
    copy_VarCoords(taC, G_zero)

   end if



   ;------normalize by global warming amounts 
   ; Set the window size for calc the mean of a period corresponding to each warming level
   window_size = 10

   ;==find the 10-yr mean corresponding to each warming level
   do ii=0,nlev-1 ;
     print("")
     print("process warming level: "+levels(ii))

     Ta_anomaly_m = Ta_anomaly(ind(models .eq. mo(mm))); select the specific model
     period_m = periods(ind(models .eq. mo(mm))); 10-yr/30-yr periods for each model
     central_year_m = central_year(ind(models .eq. mo(mm)))

     tDif = abs(levels(ii) - Ta_anomaly_m)
     ind1 = minind(tDif); closest warming amount
     yr1 = tointeger(central_year_m(ind1))
     ;matching within +-0.05 degC tolerance
     if (.not.(ismissing(ind1)) .and. tDif(ind1) .lt. 0.05) then
       print("closest matching period of "+mo(mm)+" for this level: "+period_m(ind1)+" with "+Ta_anomaly_m(ind1))
       print("at central year: "+yr1); use year index to pick the value when syr of Ta_anomaly is different from model data

      ystart = yr1 - window_size/2 +1
      yend = yr1 + window_size/2
      if (ystart.lt.syr .or. yend.gt.fyr) then
        continue
      end if
      print("start:end year for moving statistics:"+ystart+"-"+yend)

      Gday_sun_ens(mm,ii,:,:) = dim_avg_n(Gday_sun({ystart:yend},:,:), 0)
      Gday_diff_ens(mm,ii,:,:) = dim_avg_n(Gday_diff({ystart:yend},:,:), 0)
      Gday_zero_ens(mm,ii,:,:) = dim_avg_n(Gday_zero({ystart:yend},:,:), 0)

      G_sun_lev = dim_avg_n(G_sun({ystart:yend},:,:,:), 0)
      G_diff_lev = dim_avg_n(G_diff({ystart:yend},:,:,:), 0)
      G_zero_lev = dim_avg_n(G_zero({ystart:yend},:,:,:), 0)
      Z_sun_lev = dim_avg_n(Z_sun({ystart:yend},:,:,:), 0)
      Z_diff_lev = dim_avg_n(Z_diff({ystart:yend},:,:,:), 0)
      Z_zero_lev = dim_avg_n(Z_zero({ystart:yend},:,:,:), 0)

      Z_sun_lev@_FillValue = 9.96921e+36
      Z_diff_lev@_FillValue = 9.96921e+36
      Z_zero_lev@_FillValue = 9.96921e+36
      Z_sun_index = dim_pqsort_n(Z_sun_lev, 2, 0); sorting models in increasing order at each cell
      Z_diff_index = dim_pqsort_n(Z_diff_lev, 2, 0)
      Z_zero_index = dim_pqsort_n(Z_zero_lev, 2, 0)
      id_noon_sun_lev = Z_sun_index(0,:,:); index of the local solar noon
      id_noon_diff_lev = Z_diff_index(0,:,:); index of the local solar noon
      id_noon_zero_lev = Z_zero_index(0,:,:); index of the local solar noon

     G_sun_step = dim_avg_n(G_sun_lev, (/1,2/))
     print("mean G_sun at each time step (UTC time): "+G_sun_step)

     G_sun_diurnal = new((/nstep,nlat,nlon/),"float")
     G_diff_diurnal = new((/nstep,nlat,nlon/),"float")
     G_zero_diurnal = new((/nstep,nlat,nlon/),"float")
     do i=0,nlon-1
      do j=0,nlat-1
       ;daytime varies across grid cells, reorganize UTC to local time
       ;--sun scenario--
       inoon = round(id_noon_sun_lev(j,i),3)
       ;print(inoon) 
       if(any(ismissing(inoon))) then; skip missing values in oceans
          continue
       end if
       if (inoon .eq. nstep/2) then
         imorning := ispan(inoon-4,inoon-1,1)
         iaftnoon := ispan(inoon,nstep-1,1)
       else if (inoon .eq. 0) then
         imorning := ispan(inoon+4,nstep-1,1)
         iaftnoon := ispan(inoon,inoon+4-1,1)
       else if (inoon .gt. nstep/2) then
         imorning := ispan(inoon-4,inoon-1,1)
         iaftnoon := array_append_record(ispan(inoon,nstep-1,1),ispan(0,inoon-4-1,1), 0)
       else if (inoon .lt. nstep/2) then
         imorning := array_append_record(ispan(inoon+4,nstep-1,1), ispan(0,inoon-1,1), 0)
         iaftnoon := ispan(inoon,inoon+4-1,1)
       end if
       end if
       end if
       end if
       iday = array_append_record(imorning,iaftnoon, 0)
       ;print(iday)
       
       G_sun_diurnal(:,j,i) = G_sun_lev(iday,j,i) ;

       ;--shade scenario--
       inoon = round(id_noon_diff_lev(j,i),3)
       if(any(ismissing(inoon))) then; skip missing values in oceans
          continue
       end if
       if (inoon .eq. nstep/2) then
         imorning := ispan(inoon-4,inoon-1,1)
         iaftnoon := ispan(inoon,nstep-1,1)
       else if (inoon .eq. 0) then
         imorning := ispan(inoon+4,nstep-1,1)
         iaftnoon := ispan(inoon,inoon+4-1,1)
       else if (inoon .gt. nstep/2) then
         imorning := ispan(inoon-4,inoon-1,1)
         iaftnoon := array_append_record(ispan(inoon,nstep-1,1),ispan(0,inoon-4-1,1), 0)
       else if (inoon .lt. nstep/2) then
         imorning := array_append_record(ispan(inoon+4,nstep-1,1), ispan(0,inoon-1,1), 0)
         iaftnoon := ispan(inoon,inoon+4-1,1)
       end if
       end if
       end if
       end if 
       iday = array_append_record(imorning,iaftnoon, 0)
       ;print(iday)

       G_diff_diurnal(:,j,i) = G_diff_lev(iday,j,i) ; 

       ;--dark scenario--
       inoon = round(id_noon_zero_lev(j,i),3)
       if(any(ismissing(inoon))) then; skip missing values in oceans
          continue
       end if
       if (inoon .eq. nstep/2) then
         imorning := ispan(inoon-4,inoon-1,1)
         iaftnoon := ispan(inoon,nstep-1,1)
       else if (inoon .eq. 0) then
         imorning := ispan(inoon+4,nstep-1,1)
         iaftnoon := ispan(inoon,inoon+4-1,1)
       else if (inoon .gt. nstep/2) then
         imorning := ispan(inoon-4,inoon-1,1)
         iaftnoon := array_append_record(ispan(inoon,nstep-1,1),ispan(0,inoon-4-1,1), 0)
       else if (inoon .lt. nstep/2) then
         imorning := array_append_record(ispan(inoon+4,nstep-1,1), ispan(0,inoon-1,1), 0)
         iaftnoon := ispan(inoon,inoon+4-1,1)
       end if
       end if
       end if
       end if 
       iday = array_append_record(imorning,iaftnoon, 0)

       G_zero_diurnal(:,j,i) = G_zero_lev(iday,j,i) ;

     end do
    end do

      G_sun_ens(mm,ii,:,:,:) = G_sun_diurnal; 
      G_diff_ens(mm,ii,:,:,:) = G_diff_diurnal; 
      G_zero_ens(mm,ii,:,:,:) = G_zero_diurnal; 

      G_sun_step := dim_avg_n(G_sun_diurnal, (/1,2/))
      print("mean G_sun at each time step (local time): "+G_sun_step)

     end if
   end do

  end do; 12 models

   ;======calculate ensemble median and confidence interval at each grid cell at each warming level
   Gday_sun_med = dim_median_n(Gday_sun_ens, 0)
   Gday_diff_med = dim_median_n(Gday_diff_ens, 0)
   Gday_zero_med = dim_median_n(Gday_zero_ens, 0)
   G_sun_med = dim_median_n(G_sun_ens, 0)
   G_diff_med = dim_median_n(G_diff_ens, 0)
   G_zero_med = dim_median_n(G_zero_ens, 0)

   ;=========calculate confidence interval at each grid cell
   print("")
   ncrit = round(nmod/2.0, 3); round to integer 5 or 6 or 7 ;(half of models)
   print("nmod="+nmod+" | ncrit="+ncrit)

   ;show 25–75th percentiles -> 50% confidence interval (often shown in other studies, Zhang et al. 2021, Matthews et al. 2017)
   nval := dim_num_n(.not.ismissing(Gday_sun_ens),0); Count the number of True models (those provide data) at each warming level
   nval := where(nval.ge.ncrit, nval, default_fillvalue(typeof(nval))); only consider grid cells with >=6 models
   nval@_FillValue = default_fillvalue(typeof(nval))
   printMinMax(nval,1)
   pctl25_sun = round(0.25*nval,3)-1; 3rd low -> 3 models (0,1,2) at or below this value
   pctl75_sun = round(0.75*nval,3)-1; 9th high -> 9 models (0-8) at or below this value
   ;subtract 1 to get the correct model indices

   nval := dim_num_n(.not.ismissing(Gday_diff_ens),0); Count the number of True models (those provide data) at each warming level
   nval := where(nval.ge.ncrit, nval, default_fillvalue(typeof(nval))); only consider grid cells with >=6 models
   nval@_FillValue = default_fillvalue(typeof(nval))
   printMinMax(nval,1)
   pctl25_diff = round(0.25*nval,3)-1; 3rd low -> 3 models (0,1,2) at or below this value
   pctl75_diff = round(0.75*nval,3)-1; 9th high -> 9 models (0-8) at or below this value

   nval := dim_num_n(.not.ismissing(Gday_zero_ens),0); Count the number of True models (those provide data) at each warming level
   nval := where(nval.ge.ncrit, nval, default_fillvalue(typeof(nval))); only consider grid cells with >=6 models
   nval@_FillValue = default_fillvalue(typeof(nval))
   printMinMax(nval,1)
   pctl25_zero = round(0.25*nval,3)-1; 3rd low -> 3 models (0,1,2) at or below this value
   pctl75_zero = round(0.75*nval,3)-1; 9th high -> 9 models (0-8) at or below this value

   G_sun_sort = G_sun_ens
   G_sun_index = dim_pqsort_n(G_sun_sort, 2, 0); sorting models in increasing order at each cell
   G_diff_sort = G_diff_ens
   G_diff_index = dim_pqsort_n(G_diff_sort, 2, 0); sorting models in increasing order at each cell
   G_zero_sort = G_zero_ens
   G_zero_index = dim_pqsort_n(G_zero_sort, 2, 0); sorting models in increasing order at each cell

   G_sun_low = new((/nlev,nstep,nlat,nlon/), "float")
   G_sun_high = new((/nlev,nstep,nlat,nlon/), "float")
   G_diff_low = new((/nlev,nstep,nlat,nlon/), "float")
   G_diff_high = new((/nlev,nstep,nlat,nlon/), "float")
   G_zero_low = new((/nlev,nstep,nlat,nlon/), "float")
   G_zero_high = new((/nlev,nstep,nlat,nlon/), "float")
  do i=0,nlon-1
   do j=0,nlat-1
    do t=0,nstep-1
     do ii=0,nlev-1
      m_25 = pctl25_sun(ii,j,i)
      m_75 = pctl75_sun(ii,j,i)
      if(ismissing(m_25) .or. ismissing(m_75)) then;
        continue
      end if

      ip_25 = G_sun_index(m_25,ii,t,j,i); index of 25th percentile of G at each grid cell
      ip_75 = G_sun_index(m_75,ii,t,j,i); index of 75th percentile of G at each grid cell
      G_sun_low(ii,t,j,i) = G_sun_ens(ip_25,ii,t,j,i)
      G_sun_high(ii,t,j,i) = G_sun_ens(ip_75,ii,t,j,i)

      m_25 = pctl25_diff(ii,j,i)
      m_75 = pctl75_diff(ii,j,i)
      ip_25 = G_diff_index(m_25,ii,t,j,i); index of 25th percentile of G at each grid cell
      ip_75 = G_diff_index(m_75,ii,t,j,i); index of 75th percentile of G at each grid cell
      G_diff_low(ii,t,j,i) = G_diff_ens(ip_25,ii,t,j,i)
      G_diff_high(ii,t,j,i) = G_diff_ens(ip_75,ii,t,j,i)

      m_25 = pctl25_zero(ii,j,i)
      m_75 = pctl75_zero(ii,j,i)
      ip_25 = G_zero_index(m_25,ii,t,j,i); index of 25th percentile of G at each grid cell
      ip_75 = G_zero_index(m_75,ii,t,j,i); index of 75th percentile of G at each grid cell
      G_zero_low(ii,t,j,i) = G_zero_ens(ip_25,ii,t,j,i)
      G_zero_high(ii,t,j,i) = G_zero_ens(ip_75,ii,t,j,i)

     end do
    end do
   end do
  end do

  level      := levels
  level!0    = "level" ;named variables
  level@long_name = "global warming levels 1-4 degree Celsius"
  step      :=  ispan(1,nstep,1)
  step!0    = "step" ;named variables
  step@long_name = "time step of a day"
  
  G_sun_med!0="level"
  G_sun_med&level = levels 
  G_sun_med!1="step"
  G_sun_med&step = step
  G_sun_med!2="lat"
  G_sun_med&lat = lat 
  G_sun_med!3="lon"
  G_sun_med&lon = lon 
 
  copy_VarCoords(G_sun_med,G_sun_low)
  copy_VarCoords(G_sun_med,G_sun_high)
  copy_VarCoords(G_sun_med,G_diff_med)
  copy_VarCoords(G_sun_med,G_diff_low)
  copy_VarCoords(G_sun_med,G_diff_high)
  copy_VarCoords(G_sun_med,G_zero_med)
  copy_VarCoords(G_sun_med,G_zero_low)
  copy_VarCoords(G_sun_med,G_zero_high)

  Gday_sun_med!0="level"
  Gday_sun_med&level = levels 
  Gday_sun_med!1="lat"
  Gday_sun_med&lat = lat 
  Gday_sun_med!2="lon"
  Gday_sun_med&lon = lon 

  copy_VarCoords(Gday_sun_med, Gday_diff_med)
  copy_VarCoords(Gday_sun_med, Gday_zero_med)


 ;-whether to save data as ncdf 
   netCDF     = True;False
   if (netCDF) then
     dirNc  = dirPath2       ; where nc file will be saved
     filNc  = "G_3h_diel_localtime_outdoor.nc"
     pthNc  = dirNc+filNc
     system("/bin/rm -f "+pthNc)       ; remove any pre-existing file
     ncdf = addfile(pthNc ,"c")        ; open output netCDF file

     ;===================================================================
     ; create global attributes of the file (optional)
     ;===================================================================
     fAtt               = True            ; assign file attributes
     fAtt@title         = "reorganized 3-hourly data to local time according to solar zenith angle"
     if (isvar("case")) then
        fAtt@case       = case
     end if
  
     fAtt@Conventions   = "None"
     fAtt@creation_date = systemfunc ("date")
     fileattdef( ncdf, fAtt )            ; copy file attributes

     ncdf->Gday_sun_med  = Gday_sun_med
     ncdf->Gday_diff_med  = Gday_diff_med
     ncdf->Gday_zero_med  = Gday_zero_med

     ncdf->G_sun_med  = G_sun_med
     ncdf->G_sun_low  = G_sun_low
     ncdf->G_sun_high  = G_sun_high
     ncdf->G_diff_med  = G_diff_med
     ncdf->G_diff_low  = G_diff_low
     ncdf->G_diff_high  = G_diff_high
     ncdf->G_zero_med  = G_zero_med
     ncdf->G_zero_low  = G_zero_low
     ncdf->G_zero_high  = G_zero_high

   end if



 end if


end 
